Disturbance rejection model predictive control of lower limb rehabilitation exoskeleton

Nowadays, exoskeleton is broadly used in the rehabilitation training of many postoperative patients. However, the uncertainty and disturbances caused by different patients and system itself may lead to incompletely rehabilitation training as planned, or even unsafety. This paper addresses the control problem of a lower limb exoskeleton, in the spirit of the recent progress on model predictive control (MPC) and extended state observer (ESO). More precisely, our approach is based on the strategy that designing an ESO to estimate the total disturbance of the dynamics model and compensating it in the design of the MPC process. To accomplish this, we introduce the virtual control quantity to decouple the dynamics model of the system and summarize the human disturbances, unmeasured states and system non-linearity as the total disturbance of the model. By doing so, the uncertainty can be estimated by our designed ESO. Based on the moving horizontal optimization and feedback mechanism of MPC, the output prediction of the system can be more accurate since the uncertainty are effectively compensated. The virtual experiment results demonstrate that proposed controller significantly improves the control accuracy on lower limb rehabilitation exoskeleton with disturbances (improved by over 34\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\%$$\end{document}%), comparing with conventional MPC and fuzzy PID. As a result, our achievements will make contributions to better rehabilitation training for patients using rehabilitation exoskeletons.


Dynamic models
In this section, a two degrees of freedom (DoFs) lower limb exoskeleton is discussed, from which the dynamics model for the MPC controller is derived.Each limb of the exoskeleton is composed of two rotary joints and two links and is modeled as a double pendulum.The two joints are placed at the hip and knee of the wearer respectively, and the two links are fixed with thigh and shank respectively.The model of the lower limb exoskeleton is depicted in Fig. 1.It is worth noting that O point, the center of the waist, is relatively fixed in space or moving at a constant speed.This assumption is consistent with the practical application of rehabilitation exoskeletons, where the patient is relatively fixed by a bracket or sling and moves forward at a constant speed.The MPC control scheme is to make the exoskeleton taking a step following the set trajectory autonomously.To simplify the study, we only considered the swing action of one limb during a stepping cycle in sagittal plane (anterior-posterior direction).Thus, two DoFs of the motion are monitored, namely hip flexion angle θ and knee flexion angle ϕ.
The double pendulum can be considered as a potentially chaotic system, since a slight change of certain parameters in initial conditions can cause dramatic impacts on the subsequent motion.The objective of this study is to make the double links track the motion trajectory of the thigh and shank as closely as possible to ensure that the smooth and safe rehabilitation training of the patients, which means the tracking error under disturbances must be reduced to the minimum.Based on Second Lagrange equations, the dynamics of the exoskeleton can be described as, where m 1 , l 1 and J 1 represent the mass, the moment of inertia and the length of the thigh, respectively.m 2 , J 2 and l 2 represent the mass, the moment of inertia and the length of the shank, respectively.θ and θ , φ and φ are the angular velocity and angular acceleration of hip joint and knee joint, respectively.l c1 and l c2 represent the distances from the joints to the mass centers of the two legs.τ 1 and τ 2 are the torques exerted respectively on the hip joint and knee joint.τ t and τ s are the unknown torques generated by human-machine interaction.
In practice, it is reasonable to assume that the unknown control torque τ t and τ s are continuously differ- entiable.Define the state variable vector as Therefore, based on the dynamics model of the system, the exoskeleton system is presented in state-space form as, where where d(x, τ d ) are considered as the total disturbance of the system, which are caused by the human-machine interaction, nonlinear parameters deviations and so on.γ is the coefficient.They are determined as follows, It is worth mention that in this practical work, we focused on applications and used the continuous dynamical model to model the system and discretized it for more convenient use, and then, all the stability analysis is made upon the discretized version.

Model analysis
In this section, the physical parameters of exoskeleton system (i.e.m, l, J, θ , ϕ ) are utilized in the design of ESO and MPC controller, based on the MPC theory.Obviously, it can be seen from model (1) that the system is always subjected to the parametric uncertainties and disturbances (i.e. the uncertain time-varying parameters x and τ d in d(x, τ d ) ). Accordingly, we treat the total uncertainty of the model in the state-space equation as the main Each limb of the exoskeleton is composed of two rotary joints (hip and knee) and two links (thigh and shank), and modeled as a double pendulum.
disturbances, which would be estimated by the designed observer and then compensated in the MPC controller to enhance the tracking performance and robustness of the system under disturbances.In order to accomplish this task, we must first solve the problem of system decoupling.For strong coupling systems, decoupling control is an effective solution.The decoupling process 23 is realized by solving the inverse matrix using the matrix decomposition method.This method is not only simple, but also capable to deal with the time-varying matrix, so it can achieve the effect of real-time control.In addition, this method has relatively sufficient mathematical guarantee.By introducing virtual control variables: • u , the input-output relationships of link i of the exoskeleton system can be written as This means the input on link i is u D i and its output is y i = x i , so there is a single input-single output relationship between the virtual control variable u D i and the controlled output y i of each link.That is, the controlled output y i of link i and the virtual control variable u D i are completely decoupled.Therefore, the decoupling control of multivariable system can be realized by embedding n controllers in parallel between the control vector u D and the output vector y.Meanwhile, the actual control quantity u can be determined by the formula of the virtual control quantity u D Therefore, the original state-space Eq. (1) of the system can be written as, Let the sampling time of the system be T s .Notice that the sampling time for the system has to be short enough to meet the performance requirements of the exoskeleton system.Applying the forward Euler discrete method to (4), the discrete time model of the system is given by, where Ã = (I + T S • A) , B = T S • B , I T S = T S • I.
For the design of MPC controller and ESO, the following assumptions are made in this paper.

Assumption 1
The time-varying non-linear disturbances d I and d II in (1) are both bounded and satisfy where δ 1 , δ 2 , δ 3 and δ 4 are positive constant.
Assumption 2 According to the MPC theory, let the prediction horizon as N p and the control horizon of the system as N c N c ≤ N p .The state prediction equation and output prediction equation can be derived from the state-space equation based on the following assumptions, where u, d and �γ denote the increment of control law u, system disturbance d and coefficient γ respectively.

ESO design
The objective of designing ESO is to estimate the total model disturbances in real time.Accordingly, we first extend the uncertainties d I and d II as additional state variables, and let h 1 and h 2 represent the time derivative of d I and d II respectively (i.e.h I = ḋI , h II = ḋII ).The extended states are defined as 5), the extended state-space equation of the system is expressed as follows, ( www.nature.com/scientificreports/ where 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 Let z denotes the estimation of x e .From the extended system model ( 6), a linear form ESO 18 can be designed as follows, where T represents the observer gain and ω o > 0 are considered as the bandwidth of the observers.
To analyze the designed ESO 24 , the characteristic polynomial in ( 7) can be inferred to Let xe represents the estimated error (i.e.xe = x e − z ).From ( 6) and ( 7), the estimation error can be shown as follows, • Q x e , the estimation error can be rewritten as where M e = [0, 0, 1, 0, 0, 1] T and A ε is the Hurwitz matrix that can be inferred from (9).
Lemma 1 Reference 24 Assuming that H e is bounded, then the estimated states are always bounded.There exist a constant σ i > 0 , some positive integer c and a finite time T 1 > 0 , that satisfy the following equation, From Lemma 1 and (10), it can be inferred that the designed ESO in ( 7) is stable.In addition, the estimation error xe of system uncertainties can be made arbitrarily small by increasing the bandwidth ω o , according to the primary analysis 24,25 .
Therefore, the estimated values of total disturbances d I and d II can be obtained from the observed value z 3 and z 6 (i.e.dI = z 3 , dII = z 6 ).The discrete time model of the ESO is given by, where Ãe = I + T s • A e , Be = T s • B e , L = T s • L.

MPC controller design
By subtracting the states of two adjacent sampling instants, the incremental form model of the system can be derived from (5).It is expressed as, is the disturbance increment.Taking the current time step k and the current system state x(k) as the initial conditions.The state increment prediction �x(k + 1) at time step k + 1 will be represented as �x(k + 1 | k) .The predictive system output y(k + 1) at time step k + 1 will be represented as y(k + 1 | k) .Assuming that the prediction time domain is N p and the control time domain is N c .Based on Assumption 2 and (12), the state prediction equation can be conducted as follows, Vol:.( 1234567890) www.nature.com/scientificreports/Additionally, the output prediction equation can be derived based on (12) and (13).It is expressed as follows, Define Y(k) and �U(k) as follows, Therefore, the output prediction equation ( 14) can be reformed as, where the matrices of H x , H I , H d and H u are determined by the matrices of Ã, B , I T S and C in (14).The estimated error in (15) has been neglected because it can be made small enough by increasing the ESO gains.The effect of the estimation error will be analyzed in the next subsection.The objective function of MPC is designed to reflect the control and tracking performance of the controlled exoskeleton system.And the optimal control sequence can be obtained by optimizing the quadratic objective function, in order to track the desired trajectory.The objective function of the controlled exoskeleton system is designed as follows, where the first term to the right of the equal sign is the terminal constraint, and the second term is the control law increment constraint (dynamic constraint).Where R and Q in ( 16) are the diagonal weight matrices of tracking error and control increment, respectively.In addition, X r (k + 1) is an N p -dimensional vector, which represents the reference trajectory.Substituting (15) to (16), the optimal control sequence solved can be expressed as, can be calculated based on the reference trajectory and system outputs in real-time.
According to the moving horizontal optimization mechanism of MPC controller, only the first term of the optimized control sequence is input to the controlled system during each sampling period, and this optimization process will be repeated in the next sampling period.The optimal control law of MPC can be described by, where can be regarded as the predictive control gain.It is not hard to notice that the K mpc can be pre-calculated offline, which will reduce computing load of the processor.( 13) www.nature.com/scientificreports/

Stability proof
The stability of the proposed ESO-based MPC controller is analyzed in this subsection, based on the stability proof of Lyapunov theoretical framework.First, substituting y(k) = C • �x(k) + y(k − 1) to ( 17), it can be rewritten as, Substituting ( 19) to ( 12), the incremental model of the discrete system with disturbances can be written as follows, Similarly, substituting (19), 7), the incremental model of the discrete ESO can be derived as follows, Based on ( 20) and ( 21), the state increment of the closed-loop system is defined as, Consequently, the closed-loop system of ESO-based MPC can be described as, where The closed-loop system can obtain nominal asymptotic stability if the total eigenvalue of the matrix A l is within the unit circle 26 .

Remark 1
The parameters adjustment of the proposed controller and observer is discussed.For the MPC controller, the matrices R and Q in the objective function reflect the design expectation of the tracking performance and control increment.The matrix R is the weighting factor matrix of tracking error.The larger the matrix R, the closer the theoretical system output is to the reference trajectory.The matrix Q is the weighting factor matrix of control increment.The larger the matrix Q, the smaller the control quantity of the system.Meanwhile, the matrices of R and Q are also crucial to the stability of the closed-loop system.In practice, the overlarge matrix R may cause system instability due to the system state chattering.Therefore, in the process of matrix parameters adjustment, the matrix R is first adjusted based on the performance of tracking error, and then the matrix Q is adjusted to meet the limitation of the control increment and alleviate the control input chattering.
Remark 2 Another set of important parameters of the MPC controller is the prediction horizon N p and the control horizon N c .Ideally, N p should be long enough to predict the system dynamic response as long as pos- sible.However, it is unrealistic to predict the whole process due to the limitation of computational cost and time.Hence, a compromise scheme is usually adopted, which is taking the lower limit of the prediction horizon length.It is obvious that the change of u at time step k will have a direct effect on the control input u at the next time step k + 1 , and then have a long-term effect on the system response of acceleration ÿ at the time step k + 2 , on the velocity ẏ at the time step k + 3 and on the displacement y at the time step k + 4 .Therefore, we choose a minimum value of the prediction horizon N p = 5 in the implementation.Apparently, the length of the control horizon N c should be less than the prediction horizon.N c represents the number of steps which the input control quantity is seen as constant when predicting the future output of the system 27 .Similarly, to reduce the calculation time and amount of the optimization process, we choose N c = 2 in this paper.
Remark 3 ω o is considered as the bandwidth of the designed ESO.Theoretically, by increasing the bandwidth ω o , the estimation error of the total disturbance and state can be made arbitrarily small 28 .But overlarge ω o means increasing the sensitivity of the ESO to the noise, which may cause the chattering of the estimation because of the amplification of the measurement noise, etc.The adjustment process of ω 0 will be started from a small value ( 19) Vol:.( 1234567890 www.nature.com/scientificreports/and gradually increased to improve the performance of observation and the overall tracking performance until reaching a satisfied level in implementation.

Experiment and results
In this section, we explore the tracking performance of the proposed ESO-based MPC controller when applied to the task of exoskeleton motion control.Through comparing with the conventional MPC strategy, we attempt to demonstrate the advantages of the proposed controller in tracking accuracy with disturbances.

Setup and evaluation
A co-simulation model of the exoskeleton system has been built up, as shown in Fig. 2, in order to test and verify the tracking performance of the proposed controller.We firstly model the existing lower limb exoskeleton product designed by our cooperative team 29 in computer via ADAMS.Considering the calculation cost, some non-critical parts of the exoskeleton is simplified, such as screws, rubber pads, straps, etc. Besides, all joints of exoskeleton (i.e.hip, knee, ankle) are setup with CONTACT connections to simulate the real status of the system.For the power supply problem, we replace the actual driving motors (CARTRIDGE DDR Motor C043A and C053A, Kollmorgen, USA) with an ideal virtual driving force, and limit the control force increment in MPC controller based on the physical characteristics of the motor . By doing so, we avoid discuss the power and control characteristics of the motor in co-simulation.
Though previous experiments, the interaction forces between human body and the exoskeleton has been measured, as shown in Fig. 3.It is considered as external disturbances and applied to the corresponding position of the model.The reference trajectory of the human lower limbs is given in the work of other researchers 30 .Based on the composition of exoskeleton system, the nominal values of each component parameters are listed in Table 1.The control strategy, which consists of ESO, MPC controller, software interaction and data collection, is realized via MATLAB Simulink.The sampling time of the control system is set as T s = 0.01s.
We evaluate the tracking performance of our control approach and the other two control algorithms (i.e.fuzzy-PID and conventional MPC) both quantitatively and qualitatively.For quantitative description, the tracking performance is visualized in figures.Apparently, our controller show clear improves of tracking accuracy with disturbances.
For the quantitative evaluation, three performance indices, including the maximum tracking error M e , the average µ and the standard deviation σ of the tracking error, are utilized to evaluate the tracking accuracy of each method.They are defined as follows: (1) The maximum absolute value of the overall tracking errors, which is defined as, where N is the number of recorded samples.e(i) is the tracking error of each sample.
(2) The average of the tracking error, which is defined as, ( 23) www.nature.com/scientificreports/ (3) The standard deviation of the tracking error, which is defined as,

Result and discussion
Before formally tracking the trajectory of human lower limb, we first conducted a control test of hip joint angle to evaluate the speed of the proposed controller, and compared the result with fuzzy-PID and conventional MPC.
The experimental conditions are consistent with the previous description, but the reference trajectory is set as a step function from 5 degrees to 37 degrees (as shown by the solid line in Fig. 4).The performances of the three controllers are shown in Fig. 4. It can be seen that our controller is inferior to PID in terms of control speed but superior to traditional MPC, and significantly superior to both in terms of overshoot.In addition, our strategy also has better performance in terms of terminal accuracy.The control input and the overall tracking performance of the aforementioned three controllers (i.e.fuzzy-PID, conventional MPC and our proposed controller) are shown in Figs. 5 and 6, respectively.The performance indices of the three controllers are collected in Table 2.The simulation results show that our proposed controller outperforms the other two controllers on tracking performance and has stronger robustness.With the help of Moment of inertia of thigh 0.0590 Moment of inertia of shank 0.0175 the estimation by the designed ESO, our controller achieves higher tracking accuracy than conventional MPC controller, which also demonstrates the disturbances rejection ability of the designed ESO-based MPC controller.Specifically, the main reason why our controller surpasses conventional MPC is that it not only has the modelbased compensation of MPC method, but also employed a linear ESO to estimate the total disturbance (in this case, external torques and modeling uncertainties), which is integrated to the output predictive equation of the system.In other words, by building a novel output prediction equation with ESO-estimated total disturbance, the response of the controlled exoskeleton system will be predicted more accurately, which enables the optimization solution of the control quantity to be obtained more precisely by solving the novel objective function containing disturbances.On the other hand, comparing Figs. 5 and 6, although the performance of conventional MPC  controller is not as good as ours, it is better than fuzzy-PID, which indicates that the traditional MPC method still has acceptable robustness.Besides, the control inputs of the proposed controllers are bounded, as shown in the Fig. 7.
As shown in Table 2, our controller shows better tracking accuracy in all indices.By introducing linear ESO to the MPC, the mean tracking error of the proposed controller is reduced down to 0.93 degree (hip flexion angle) and 1.07 degree (knee flexion angle), while the mean tracking error of MPC controller is about 1.45 degree (hip flexion angle) and 1.61 degree (knee flexion angle).Compared with conventional MPC, the tracking performance of our controller is improved by about 34-35% .Furthermore, our controller has greater advantage compared with fuzzy-PID, the tracking performance is improved by about 35-39%.

Conclusions
In this paper, we have proposed a disturbance rejection model predictive control strategy based on ESO for a lower limb rehabilitation exoskeleton.Our key is to integrate the estimation of system total disturbance to the quadratic objective function, which is used for solving the optimal control law of MPC.To do so, we first introduce virtual control variables to decouple the exoskeleton system, and then construct an ESO based on the expanded decoupled system, that can estimate the model uncertainties and the disturbance from the wearer as a total disturbance.Based on the stability theory of Lyapunov framework, the stability of the closed-loop system is ensured, which indicates that the proposed ESO-based MPC controller has prescribed tracking performance under the condition of disturbances and uncertainties.Through comparative virtual experimental results, we have demonstrated that our control scheme for exoskeleton outperforms two existing techniques (fuzzy-PID and conventional MPC) on the effectiveness and priority, in the strong disturbance case.As a result, our achievements will contribute to better rehabilitation training of patients using lower limb rehabilitation exoskeleton.

Limitations and future works
In this paper, a certain simplification has been made to the physical model of the lower limb rehabilitation exoskeleton, and the control strategies have been proposed just for the gait in sagittal plane, to simplify the study.In addition, we regard the human-machine interaction forces and system parameter uncertainty as total disturbances to be observed and compensated, which poses poses a challenge to the performance of ESO.In future works, further refinement of the physical model of the exoskeleton can be carried out, and the problem of gait walking in three-dimensional space can be considered.In addition, the design of ESO can further observe interaction forces and uncertainties separately, depending on how the system is modeled.

Figure 1 .
Figure 1.Illustrations of the lower limb exoskeleton (left) and the dynamics model in sagittal plane (right).Each limb of the exoskeleton is composed of two rotary joints (hip and knee) and two links (thigh and shank), and modeled as a double pendulum. https://doi.org/10.1038/s41598-023-46885-4

Figure 2 .
Figure 2. Control flow diagram of the co-simulation model.The deviation between the ESO estimated position and the reference position, combining with the total disturbance also obtained by ESO estimation, is substituted into the MPC controller to obtain the optimal control torque.The control quantity is input into the exoskeleton dynamics model to obtain the actual position, which is feed back to the controller to form a closed loop control system.

Figure 3 .
Figure 3.The interaction forces between human body and the exoskeleton.The upper sub-figure is the interaction force between thigh and the exoskeleton.The lower sub-figure is the interaction force between shank and the exoskeleton.

Figure 4 .
Figure 4.The tracking performance of three control scheme.The solid line is the reference trajectory we set as a step function.The dotted line represents the PID control result, the dashed line represents the MPC control result, and the dash-doted line represents our strategy control result.

Figure 5 .
Figure 5.The tracking performance of three schemes for hip joint control.The solid line is the reference trajectory, which refers to the hip flexion angle of the human hip joint during one gait cycle.The dotted line represents the fuzzy-PID control result, the dashed line represents the conventional MPC control result, and the dash-doted line represents our strategy control result.

Figure 6 .
Figure 6.The tracking performance of three schemes for knee joint control.The solid line is the reference trajectory, which refers to the knee flexion angle of the human knee joint during one gait cycle.The dotted line represents the fuzzy-PID control result, the dashed line represents the conventional MPC control result, and the dash-doted line represents our strategy control result.

Figure 7 .
Figure 7.The control inputs of our control scheme.The upper sub-figure is the control torque on hip joint and the lower sub-figure is the control torque on knee joint.

Table 1 .
Parameters of the exoskeleton system.

Table 2 .
Performance indices of the aforementioned three control schemes.